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ABSTRACT 

A number of globular clusters appear to have undergone core collapse, in the sense 
that their predicted collapse time is much shorter than their current age. Simulations 
using gas models and Fokker-Planck approximation have shown that the central density 
of a globular cluster after the collapse undergoes nonlinear oscillation with large am- 
plitude (gravothermal oscillation). However, whether such an oscillation actually takes 
place in a real iV-body system has remained unsolved, because an iV-body simulation 
with a sufficiently high resolution would have required the computing resource of the 
order of several Gflops-years. In the present paper, we report the result of such a sim- 
ulation, performed on a dedicated special-purpose computer GRAPE-4. We simulated 
the evolution of isolated point-mass systems with up to 32,768 particles. The largest 
number of particles reported previously is 10,000. We confirmed that gravothermal os- 
cillation takes place in an iV-body system. The expansion phase shows all signatures 
that are considered as the evidences of the gravothermal nature of the oscillation. At 
the maximum expansion, the core radius is ~ 1% of the half-mass radius for the run 
with 32,768 particles. The maximum core size r c depends on N, as < r c >oc TV -1 / 3 . 



1. Introduction 

A number of globul ar clusters appear to have un- 
dergone core co llapse ( Cohn and Hut 1984 , Hut and 
Djogovski 1992j ). Djorgovski and King (1986) showed 
that 15% of galactic globular clusters have unresolved 
density cusps. Recent observations with HST have 
demonstrated that some of these clusters have the 
core sizes smaller than 0.03 arcsec, which i s as far as 
one can go with present technique (e.g., sosin and 
King 1995| , |Yanny et al. 1994[ )- 



The theoretical study of the evolution after the 
core collapse was pioneered by Henon (1975), who 
incorporated the energy production by binaries into 
his Monte-Carlo calculation. He found that the 
whole cluster expanded homologously in the thermal 
timescale. Similar results were obtained by gaseous 
models and Fokker-Planck calculations. 

Sugimoto and Bettwieser (1983, Bettwieser and 
Sugimoto 1984) found that the post-collapse expan- 
sion is unstable if the energy production is inefficient. 
In the unstable case, the central density showed an 
oscillation with a large amplitude. They called this 
oscillation gravothermal oscillation. They modelled 
the post-collapse evolution of globular clusters using 
a conducting gas sphere with artificial energy pro- 
duction. The efficiency of the energy production is 
related to the total number of particles N, and for 
large N efficiency is small. Thus, their result implies 
an iV-body system with number of particles larger 
than a critical number should show the gravothermal 
oscillation. 

Since several other calculations based on similar 
models did not reproduce the result of Sugimoto and 
Bettwieser (1983), the validity of their result had 
been controversial for a few years after their first 
paper. However, calculations with improved accu- 
racy confirmed that the oscillation takes place in both 
the g aseous models a n d Fokker-Planck (FP) calcula- 
tions( Goodman 1987 , Hcggic and Ramamani 1989] , 
Cohn et al. 1989|) . 



What have been found in gas models and FP cal- 



culations are summarized as follows. There 
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as ionows. mere is a 
mologously expanding solution for the posr-collapse 



evolution of the cluster ( Goodman 1984 ). In this so- 
lution, the size of the core is proportional to N^ 1 ^ 3 . 
Thus, the core is smaller for larger number of parti- 
cles. If the core is too small, this homologous expan- 
sion becomes unstable in a way similar to the way in 
which an isothermal sphere is unstable, and the core 



starts to show an oscillating behavior. The critical 
number of particles is around 7,000 (Goodman 1987). 
If the total number of particle is slightly larger than 
this critical value, the oscillation is regular. However, 
at TV ~ 10,000, "period doubling" takes place and 
the density shows oscillation with double peaks. For 
N > 50, 000, oscillation becomes apparently chaotic 



(Heggie and Ramamani 1989, Breeden et al. 1994). 

Both the contraction and the expansion are driven 
by the gravothermal instability, and the energy pro- 
duction by binaries acts as the trigger to start the 
expansion. During the expansion, the thermal energy 
flows in from the halo to the core, and this heat flux is 
supported by so-called "temperature inversion" . This 
temperature inversion itself is formed through the ex- 
pansion of the core driven by the binary heating in 
the following way. 

Binaries deposit the energy to the core through 
"indirect heating" ( Hut 1985 ). When a binary hard- 
ens through interaction with a third star, that star is 
likely to be kicked out of the core, because the recoil 
velocity is much larger than the velocity dispersion of 
the core. As a result, the binding energy of the core 
becomes smaller. The core expands to reestablish the 
dynamical equilibrium. This expansion causes the ve- 
locity dispersion of the core to decrease. Thus the 
central temperature becomes lower than the temper- 
ature of the region just outside the core. 

Once the central temperature becomes lower than 
the temperature around the core, the heat starts to 
flow inward. This inward heat flux let the core expand 
further, and the temperature of the core continues to 
decrease. This self-supported expansion can continue 
as far as the temperature gradient outside the core 
is small. In practice, however, temperature gradient 
becomes larger as the core size becomes larger, and 
the temperature inversion vanishes at a certain core 
size. After that, the core starts to recollapse. 

The critical number of particles for multi-component 
systems is larger than that for single-component sys- 
tems, and for a "realistic" mass function, the critical 



number of particle can be as large as 3 x 10 (Murphy 



et al. 1990) 



The question whether such gravothermal oscilla- 
tion would takes place in real globular clusters or 
in A^-body systems has not been settled yet. Both 
gaseous models and FP models have many simplify- 
ing assumptions that might make the evolution com- 
pletely different from that of a real cluster. For ex- 
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ample, both assume spherical symmetry, while in TV- 
body simulations cores are known to wander around 
flMakino and Sugimoto 1987| , |Heggie et al. 1994[ ). 
Both assume that the energy production by binaries is 
smooth, while in TV-body system binaries are formed 
stochastically and emit energy intermittently. 

An TV-body simulation with a sufficiently large 
number of particles has been impossible, simply be- 
cause the requirement for computer power has been 
excessive. Makino et al. (1986) performed 100-body 
simulations of an equal-mass cluster, using both soft- 
ened and unsoftened potentials. They used a soft- 
ened potential in order to reduce the energy produc- 
tion from binaries. They found oscillatory behavior 
for both cases, but with larger amplitude for softened 
potential. They argued that the fact that the ampli- 
tude of oscillation is larger for reduced energy pro- 
duction implies that the observed oscillation is driven 
by the gravothermal instability and not purely by the 
heating by binar y. McMillan and Lightman (1984, 
McMillan 1986a ) developed hybrid calculation code 
in which the central region of a cluster is treated as 
an JV-body system while the outer region is treated 
in FP formalism. They observed oscillatory behav- 
ior but concluded that it was not gravothermal. In- 
agaki (1986) performed 1000- and 3000-body simula- 
tions but found no sign of oscillation. He used stan- 
dard Salpeter mass function as the distribution of the 
masses of particles. Makino (1989) performed sim- 
ulation of a 3000-body equal-mass system and saw 
some signs of oscillation, including the temperature 
inversion. Spurzem and Aarseth (1996) performed 
10,000 body simulation and saw oscillatory behav- 
ior, but they did not see any clear evidence of the 
gravothermal oscillation. TV-body simulations in the 
last decade can be summarized as follows: The cen- 
tral density showed some oscillatory behavior, and 
its amplitude is seemingly larger for larger number 
of particles. However, each expansion phase is as- 
sociated with large energy input from one or a few 
binaries. As a result, whether the expansion is driven 
by instability or by energy input has been unclear. In 
addition, in all but one simulations, the number of 
particles was smaller than the critical value of 7,000. 
The 10,000 body simulation by Spurzem and Aarseth 
(1996) was not long enough to draw clear conclusion. 
Therefore it was not clear whether oscillatory behav- 
ior is caused by gravothermal instability or not. 



if one constructs an TV-body system that has an ini- 
tial temperature inversion. They constructed an TV- 
body system which mimics the density and temper- 
ature structures of the expansion phase of the gas 
model calculation by Heggie and Ramamani (1989), 
and follow the evolution of the TV-body system. They 
found that the core expands in the thermal timescale. 
Thus, they at least proved that the core of an TV-body 
system can expand gravothermally. 

In the present paper, we describe the result of 
TV-body simulations with number of particles 2,048- 
32,768. This is the first calculation with the number 
of particles well beyond the critical number of 7,000 
which covers sufficiently long time after the collapse 
to determine the nature of the post-collapse evolu- 
t ion. All simula tions were performed on GRAPE-4 
( Taiji et al. 1996 ), a special-purpose computer for col- 
lisional TV-body simulations. Our main results are the 
following. First, for large TV, the core density and 
core mass exhibited oscillations with large amplitude. 
The core mass at the maximum expansion is in good 
agreement with FP or gas model results (1% of the 
total mass for 16k- and 32k-particle runs). Second, 
we confirmed that the observed oscillation is driven 
by gravothermal instability. Several long expanding 
periods without any energy input were observed. The 
temperature inversion was visible in such expansion 
phases. The behavior of the core density is strikingly 
similar to the result of FP calculations with stochas- 



Heggie (1989, see also Heggie et al. 1994 ) tried to 
see whether gravothermal expansion would take place 



tic heat source (Takahashi and Inagaki 1991), which 
strongly suggests that the mechanism is the same. 
The trajectory in the central density-central temper- 
ature plane clearly indicates that the inward heat flux 
supports the expansion. 

The structure of this paper is the following. In 
section 2, we describe the initial model, numerical 
method and the computer used. In section 3 we 
present the result. Section 4 is for discussion. 

2. Model and Numerical method 

2.1. Initial models and the system of units 

We followed the evolution of isolated systems of 
point-mass particles. For all calculations, we used 
random realizations of the Plummcr model as the ini- 
tial condition. The system of units we adopted is the 
standard (Heggie) units ( Heggie and Mathicu 1986| ), 
in which G = 1, M = 1, and E = -1/4, where G 
is the gravitational constant, M is the total mass of 
the cluster, and E is the total energy of the cluster. 
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All particles have the same mass m = l/N, where 
N is the total number of particles. In this unit, the 
half-mass crossing time the is 2y/2. 

2.2. Numerical method 

For all calculations, we used NBODY4 (Aarseth 
1985, 1996), modified for GRAPE-4 (Taiji et al. 1996). 
The numerical integration scheme adopted in NBODY4 
is the 4th order Hcrmite scheme ( Makino 1991a| . 
Makino and Aarseth 1992) with the individual (hier- 
archical) timcstep algorithm ( McMillan 1986b| , Makino 
1991b| ) to use the GRAPE hardware or parallel/ vector 
computers efficiently. Close two-body encounters 
and stable binaries are handled by KS regularization 
(Kustaanhcimo and Stiefel 1961:). Special treatment 



for compact few-body subsystems is also possible. 

The 4th-order Hermite scheme has several advan- 
tages over traditional ABM type schemes. It is a self- 
starting scheme and therefore much easier to imple- 
ment than the multistep ABM scheme of the same 
order. Its local truncation error has much smaller co- 
efficient, which allows a larger timestep for the same 
accuracy. On general-purpose computers, this ad- 
vantage is partly cancelled by the additional CPU 
time to calculate the first time derivative, but on 
GRAPE-type machines this is an significant advan- 
tage ( [Makino and Aarseth 1992| , [Makino et al. 1993[ ). 

The output are taken at intervals of a fixed time, 
which is a certain fraction of the crossing time. We 
recorded central density, core radius, number of par- 
ticles in the core, radii of Lagrangian shells, velocity 
dispersion within Lagrangian shells, binaries and their 
binding energies. The core parameters are calculated 
following Casertano and Hut (1985), with modifica- 
tions described in McMillan et al. (1990). 

The accuracy parameter of the time integration is 
adjusted so that the energy error between two outputs 
is smaller than a certain prescribed value. The value 
we used is 1 x 10 -5 ~ 1 x 10~ 6 depending on N. 
We required higher accuracy for larger A, since the 
duration of the simulation is longer. When the energy 
error is very large, the program automatically reads 
the output at the previous checkpoint and restarts 
with a reduced accuracy parameter. 

The A-body simulation of point-mass systems with 
a large number of particles poses several technical 
problems beside the necessary computing power. The 
difficulty comes from the wide range of the timescale. 
The critical separation of two particles at which we 



apply the KS regularization is about l/N. The 
timcstep for particles involved is of the order of 
l/(100A r15 ). On the other hand, the timescale of the 
evolution of the system is of the order of N . Thus, the 
timesteps of particles that are integrated in KS for- 
malism is 1/1007V 2 - 5 of the system time, which, for 
N ~ 10 5 , about as small as the limit of double pre- 
cision numbers. Note that this estimate is valid for 
a soft binary. A hard binary would require timesteps 
smaller by several orders of magnitude than described 
above. Of course, a clever treatment can overcome 
this kind of difficulties, but the amount of program- 
ming work and the complexity of the resulting code 
would be considerable. The simplest solution would 
be to use the longer number format (e.g., the quadru- 
ple precision). This would have been a viable solu- 
tion in '70s, when mainframe machines could handle 
quadruple precision arithmetic reasonably fast. Un- 
fortunately, most of present RISC computers are ex- 
tremely slow at quadruple precision arithmetic. 

2.3. Hardware and calculation cost 



For all calculations, we used GRAPE-4 (Taiji e\ 
al. 1996| ). The GRAPE-4 is a special -purpose com- 



puter designed to accelerate the A-body simulation 
using the Hermite integrator and hierarchical timestep 
algorithm. The total system consists of 1692 pipeline 
processor chips and has the theoretical peak speed of 
1.08 Tflops. The simulations reported in the present 
paper were performed while the assemblage and test- 
ing of the GRAPE-4 system were underway. Thus 
the number of processors varies during the calcula- 
tion. For most of the 32k particle run, we used one 
quarter of the system with 423 processor chips, which 
has the peak speed of 270 Gflops. Calculations with 
smaller N were performed while the available number 
of processors was smaller. 

The actual performance depends on many fac- 
tors, but most strongly on the number of particles. 
With the present host computer (DEC 3000/900), av- 
erage speed we got for 32k-particle run was about 
50 Gflops. This 32k-particle run took about three 
CPU months. In comparison, the lOk-particle run 
reported by Spurzem and Aarseth (1996) took two 
CPU months on a Cray YMP. Since the calculation 
cost of globular cluster simulation is proportional to 
TV 3 - 3 , one quarter of GRAPE-4 is effectively 50 times 
faster than a Cray YMP. 

At present, a 32k-particle system is about the 
largest system we can try for equal-mass case. Since 
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the calculation speed of GRAPE-4 is still limited by 
the speed of the host computer, the actual calculation 
speed is roughly proportional to A for TV < 2 x 10 5 . 
Thus, the CPU time for globular cluster simulation 
on GRAPE-4 is proportional to A 2 - 3 . A 64k-particle 
run would take about one year. 

The above estimate is for an equal-mass isolated 
cluster of point-mass particles. The computational 
cost of more realistic calculations depends on many 
factors. Among them, the mass spectrum and the 



presence of primordial binaries are most important, 



and have adverse effects. 



Systems with a mass spectrum evolves faster than 



equal-mass systems (see, e.g., Inagaki and Wiyanto 



19^, |Murphy and Cohn 19881 ). Therefore it is pos- 
sible to handle a larger system. For example, if the 
evolution of the system is p times faster, GRAPE-4 
can handle the number of particles ^/p times larger. 
Thus, if the collapse is 4 times faster, GRAPE-4 can 
finish a 32k-particle run in 3 weeks or 64k-particle run 
in 3 months. 

The presence of the primordial binaries would in- 
crease the CPU time. At present, however, actual 
cost is difficult to estimate because the calculation 
cost depends strongly on the core size, of which we 
have rather little knowledge (see section 4) . It should 
be noted that the number of floating point opera- 
tions used to follow the life of a binary depends only 
weakly on the total number of particles in the system 
dMakino and Hut 1990[ ). Therefore, while the CPU 
time for the orbit integration of single stars would in- 
crease as A 2 - 3 , the cost of binaries would increases 
only as A. It is likely that for large N the cost of 
binaries is relatively small, even for simulations with 
primordial binaries. 

3. Result 



3.1. The core oscillation 

Figure 1 shows the time evolution of the central 
density for all runs. The time is scaled so that ther- 
mal timescale is the same for all runs. The scaling fac- 
tor is t r (W00)/t r (N) = 212.75 log(0.11A)/A (Giersz 
and Heggie 1994). The core density shows oscillations 
with a large amplitude in calculations with large A 
(>16k). No matter what is the real nature of this 
oscillation, it is at least clear that the core density of 
the A-body system shows oscillation with the ampli- 
tude comparable to that observed in gas models or 
FP calculations. 



For runs with small A (2k and 4k), there are some 
oscillation-like features but they are hardly distin- 
guishable from fluctuations. This result is similar to 
that of previous studies (e.g. Makino 1989). For large 
values of A, however, the oscillation with large am- 
plitude is clearly visible. 

Note that in figure 1 there is no clear transition 
from stable expansion to oscillation or from regu- 
lar oscillation to chaotic oscillation, which were ob- 



served in gas and FP models (Heggie and Ramamani 
1989|, [Brecdcn et al. 1994|). The reason is that bi 



naries emit energy intermittently (Makino and Sug- 



imoto 1987). Takahashi and Inagaki (1991) incor- 



porated this stochastic nature of the energy source 
to their FP model and found that the core density 
shows chaotic oscillatory behavior even if the energy 
production rate is larger than the critical value for 
transition between stable expansion and regular os- 
cillation. They also found that the amplitude of the 
oscillation is smaller for larger energy input (smaller 
A), which is consistent with the present result. In 
fact, it would be difficult to distinguish between our 
A-body calculation results from their stochastic FP 
results. Only visible difference is that their result is 
smoother while the central density is low. This is be- 
cause the FP result is perfectly smooth as far as there 
is no binary heating. 

3.2. Analysis of the 32k run 

In the previous section, we overview the post- 
collapse evolution of A-body point-mass systems with 
2k-32k particles. In this section, we take a closer look 
of the 32k-particle simulation to see whether we can 
find direct signature of the gravothermal expansion. 
It is generally believed that a long expansion phase 
without significant energy input and a temperature 
inversion during the long expansion phase are the 
most direct signatures of the gravothermal expansion 



(Bcttwieser and Sugimoto 1984, McMillan and Engle 
1996| ). In this section we investigate both of them. 



Figure 2 shows an enlarged view of the time vari- 
ation of the core radius as compared with the sum of 
binding energies of all binaries. It is clear that most 
of the energy is generated when the core is very small. 

The continued expansion without energy input 
from binaries has been considered to be the most 
direct signature of the gravothermal oscillation. In 
figure 2, we clearly see three such expansions, for 
t = 6700-6740, t = 6860-6920, and* = 7220-7280. 
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Fig. 1. — The logarithm of the central density plotted 
as a function of the scaled iV-body time. Curves for 
different values of N are vertically shifted by 3 units. 
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Fig. 2. — The core radius (top) and the binding en- 
ergy of binaries (bottom) as a function of time for the 
32k-particle run. 



All these expansions continue for more than 10 half- 
mass crossing times, which is hundreds of the core 
relaxation time. These expansions cannot be driven 
simply by binary heating. If the expansion were 
driven only by binary heating, it could not continue 
without energy input for the timescale more than 100 
times longer than the core relaxation time. 

Figure 3 shows temperature profiles for the con- 
tracting and expanding phases. Near the end of the 
expanding phase the temperature inversion of the or- 
der of 5% is clearly visible. Since these profiles are 
time-averaged over 10 time units (80 snapshots), the 
actual inversion might be somewhat stronger. Note 
that the temperature inversion is visible only near the 
end of long e xpansion phases also in gas model and FP 
calculati ons ( Bettwieser and Sugimoto 1984 , Cohn e; 



al. 1989| ). 



Figure 4 shows the relation between the central 
density and the central velocity dispersion. The tra- 
jectory shows striking resemblance to what obtained 
by gas model calculation (Goodman 1987). The fact 
that the trajectory shows clockwise rotations means 
that this is a refrigeration cycle, in which the central 
region absorbs heat when the temperature is low, and 
release heat when the temperature is high (Bettwieser 
and Sugimoto 1984, Bettwieser 1985). In particular, 
the later phase of the large expansions (indicated by 
the arrow marked "B") is nearly isothermal. Since 
the binding energy of binaries is unchanged during 
this phase, this nearly isothermal expansion is driven 
by the heat supplied from outside the core. In other 
words, the expansion is gravothermal. 



1.4 
1.3 
1.2 
1.1 
1 



Fig. 3. — Velocity dispersion profiles for (a) contract- 
ing and (b) expanding phases. Each profile is ob- 
tained by time averaging over 80 snapshots (10 time 
units). The time interval between curves is 5 time 
units. 
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3.3. The core size 
tion 



result and interpreta- 



Figurc 5 shows the evolution of the number of par- 
ticles in the core. For all runs, the number of parti- 
cles at the maximum contraction is of the order of 10, 
while that at the maximum expansion is ~ 1% of the 
total number of particles. This result is again in good 
agreement with gas models and FP calculations (e.g., 



Heggie and Ramamani 1989 Brecdcn et al. 1994 ). 

Figure 6 shows the fraction of time for which the 
number of particles in the core is smaller than the 
value N c as a function of N c /N, for the post-collapse 
phase. In other words, this figure shows the cumu- 
lative distribution of the core mass. For N > 8192, 
the median core mass is around 0.5-0.6 %. This cor- 
responds to r c /rh ~ 0.01. For N = 32, 768, the core 
mass is somewhat smaller than that for 16k or 8k 
runs. 

Whether the core size in the gravothermal oscilla- 
tion phase depends on the total number of particles or 
not is not well understood. Bettwieser and Sugimoto 
(1984) argued that clusters spend most of the time 
in the most expanded state, and that the structure 
of the cluster at the maximum expansion does not 
depend on the total number of particles. However, 
it is clear from figure 6 that the cumulative time is 
roughly proportional to the core radius, and that the 
coefficient is larger for larger N . In other words, the 
size of the core depends on N. 

Simple theoretical argument suggests that the core 
size must depend on N . The time-averaged energy 
production in the core would be the same as that for 



the stable exp ansion ( Goodman 1987 , Heggie and Ra- 
ma nam 1989| ) . This requirement pose a constraint on 



the time-averaged core size. The energy production 
rate is expressed as 



dE w 2 
-A " McP " 



(1) 



where M c and p c are the core mass and core density. 
Here we ignored the dependence on the velocity dis- 
persion, since the inner part is almost isothermal. For 
M c and and the core radius, r c , we have the relation 
M c cx r c , again because the inner part of the cluster 
is almost isothermal. Thus we have 



dE 
~dt 



cx m: 



(2) 



The time-averaged energy production rate must be 
equal to that of the homologous expansion solution 



o 




Fig. 4. — The change of the central density and the 
central velocity dispersion. Each data points is time- 
averaged value over 80 snapshots. Arrows indicate 
the direction of evolution. 
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Fig. 5. — The number of particles in the core as the 
function of the scaled time, for simulations with 2k, 
4k, 16k, and 32k particles 
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obtained by (Goodman 1984). The core size for 
the homologous expansion is N~ 2 ' 3 (Goodman 1984, 
Giersz and Heggie 1993). Thus we obtain 

< M~ 3 >" 1/3 oc N~ 2/3 , (3) 

where < x > means the time averaging. 

Figure 7 shows the arithmetic mean and the max- 
imum of the core mass after the collapse as functions 
of N. We also plot the quantity M e =< M~ 3 >~ 1/3 , 
since we have the theoretical prediction only for this 
quantity. For maximum core mass, the FP result by 
Breeden et al. (1994) are plotted for N — 5 x 10 4 and 
10 6 . These values are read by eye from their figure 5. 
The quantity M e indeed shows reasonable agreement 
with the theoretical prediction of equation (j|). 

However, the average and maximum core radii 
shows the dependence noticeably weaker than the the- 
oretical prediction. It should be noted that our N- 
body results for the maximum core mass and the FP 
results of Breeden et al. (1994) are on one line ex- 
pressed as M c ex N^ 1 / 3 . From figure 7, the core mass 
at the maximum expansion is approximated as 

/ N V 1/3 
M c , max ~ 0.006 ( — J M. (4) 

The slope of the average core size seems to be some- 
what steeper. 

Qualitatively, it is natural that the ratio between 
the maximum core mass and M e is larger for larger 
N, since the ratio between the minimum core mass 
and M e is larger for larger N. The minimum core 
mass is 0(1 /N), since the number of particles in the 
core at the core bounce is always around 10-30. In 
other words, the energy production rate at the core 
bounce is larger than the time-averaged rate by a fac- 
tor proportional to N. Thus, roughly speaking, the 
core mass has to stay at the value larger than the 
value at the core bounce for most of time. The typical 
relaxation timescale of the core not at the maximum 
contraction is, therefore, N times longer than that at 
the maximum contraction. Since the relaxation time 
is proportional to M 2 , this implies that the typical 
core mass is N~ x / 2 . This is in fact a reasonable fit 
for the average core mass in figure 7. 

One interesting question is whether it is possible 
to distinguish the core in the gravothermal oscillation 
and the core dominated by primordial binaries. The 
theoretical prediction by Goodman and Hut (1989) 
gives r c /rh ~ 0.02. The Fokker-Planck calculation by 




N c /N 

Fig. 6. — Fraction of time for which the number of 
particles in the core is smaller than N c , as a function 
of N c /N. 




Fig. 7. — Core mass plotted as a function of the to- 
tal number of particles N. Filled squares are M e (see 
text). Filled crosses are arithmetic means. Filled cir- 
cles are maximum values from present iV-body sim- 
ulations. Open circles are maximum values from FP 
calculations of Breeden et al. (1994). Straight lines 
indicates M c cx N~ x / 3 (top) and M c cx iV" 1 / 3 (bot- 
tom), respectively. 
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Gao et al. (1991) gave similar result. On the other 
hand, TV-body simulations by McMillan et al. (1990) 
gave r c /r h ~ 0.2. McMillan et al. (1990) used 1136 
particles. If their TV-body results can be extrapolated 
to larger values of TV, cores with primordial binaries 
and cores in gravothermal oscillation would be clearly 
distinguishable. On the other hand, FP result by Gao 
et al. (1991) indicates that the typical core size at 
the "primordial binary burning" phase and that at 
the gravothermal oscillation phase would not be much 
different. To obtain a definitive answer, we need to 
perform simulations of clusters with primordial bina- 
ries using the number of particles larger than that 
employed by McMillan et al. (1990). 

3.4. The core size — comparison with obser- 
vations 

Djorgovski and King (1986) showed that 15% of 
galactic globular clusters have unresolved density cusp. 
They classified these clusters as "post core collapse" 



PCC). Recent ground-based observations (Lugger et 
al. 1995|) and HST resu l ts (|5osin and King 1995j , 
Guhathakurta et al. 1994, Yanny et al. 1994| ) demon- 
strated that some of the PCC clusters have the cusps 
unresolved even with HST (NGC6624, M15). For 
these clusters, the surface density profile shows the 
power law down to 0.3 arcsec. Other clusters (M30), 
though the slope seems to levels off toward the center, 
turned out to be difficult to determine the slope ac- 
curately because there are too few stars. In the case 
of M30, the core radius can be anywhere smaller than 
1.5 arcsec. 

Many researchers claimed that there is a good 
agreement between the theory and observation of the 
core radius. For example, Djorgovski and Meylan 
(1994) found that no observed cluster has the core 
radius smaller than 1% of the half-light radius, using 
the data compiled by Trager et al. (1993), and argued 
this to be in good agreement with the model that as- 
sumes the binary-dominated core for PCC clusters 
(Vesperini and Chernoff 1994). 

This apparently good agreement between the "ob- 
servation" and "theory" is actually due to the fact 
that both of them overestimated the core size. On 
the observational side, High-resolution observations 
have shown that the real core size of the PCC clus- 
ters is significantly smaller than what is assumed in 
Djorgovski and Meylan (1994). Trager et al. (1993) 
adopted the HWHM as the core size of PCC clusters. 
Since the HWHM has no relation to the real core 



radius of the PCC clusters, it is quite natural that 
observations with higher resolution obtained smaller 
core radii. On the theoretical side, our present re- 
sult suggest that the typical size of the core in the 
gravothermal oscillation phase is significantly smaller 
than the previous claims. 

Our result suggests that the typical size of the core 
in the gravothermal oscillation phase is smaller than 
previous estimates because of the following two rea- 
sons. First, the time-averaged core size is about a half 
of the maximum size. Second, the core size depends 
on TV as TV _1 / Q , with index alpha 2^3. 

For NGC6624 and M15, core sizes seem to be too 
small for a stably expanding cluster, no matter what 
is the heat source. So it is highly likely that they are 
undergoing the gravothermal oscillation. 

4. Discussion 

We performed direct TV-body simulation of the 
post-collapse evolution of globular clusters. We con- 
firmed that the gravothermal oscillation actually takes 
place in a point-mass TV-body system. 

Whether real globular clusters undergo gravother- 
mal oscillation or not is a question which requires fur- 
ther research. If clusters contain many primordial bi- 
naries, even after the core collapse they might still 
be burning the primordial binaries. In addition, the 
effect of two-body binaries to the evolution of the clus- 
ter is still unclear. In fact, the cross section for the 
binary formation by tidal capture is not fully under- 
stood yet (Mardling 1995a, 1995b). 

The most straightforward way to study the effect 
of primordial binaries or two-body capture is the di- 
rect TV-body simulation. In principle, we can put pri- 
mordial binaries and their evolution into FP calcula- 
tion as Gao et al. (1991) did. However, the standard 
one-dimensional FP calculation, which assumes the 
isotropic distribution, is not appropriate to follow the 
binary population, since most binaries are formed in 
the core and their orbits are nearly radial. Thus, we 
have to solve the FP equation at least in three dimen- 
sions (E, J and the binary binding energy E^). It 
would require prohibitively large computer power. In 
addition, reliable two-dimensional FP calculation has 
become possible only very recently (Takahashi 1995). 
One could also use the Monte-Carlo approach, but 
its result must be compared with TV-body simulation 
anyway. 

We demonstrated that the TV-body simulation with 
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the number of particles close to real globular clus- 
ters is now possible, thanks to the extremely powerful 
special-purpose computer and its full-time availabil- 
ity. We are now able to use the direct TV-body simula- 
tion to study various aspects of evolution of globular 
clusters. 

Our 32k-particle calculation took about three months 
of CPU time on 1/4 of GRAPE-4. If we tried to do 
a similar calculation on a Cray T90 vector supercom- 
puter, it would have taken several years of CPU time. 
It is simply impossible to do on present-day supercom- 
puters. If we want to finish the calculation in, say, one 
CPU month, we need a computer 50-100 times faster 
than a Cray T90, which will be available in 10 years 
from now. 

The CPU time of three months is still very long. 
However, for many simulations, we do not need as 
many as 32k particles. A 16k-particle calculation was 
finished in 2-3 weeks, on 1/8 of GRAPE-4. Thus to 
run many simulations of 16k particle system is now 
practical. 

If we can continue the development of the special- 
purpose computer, we will have a system 100-1000 
times faster than present GRAPE-4 in the next five 
years. Such a system will make it possible to run 
50-100k particle simulation routinely, while 500k-lM 
particle simulation will still take months. 

I thank Sverre Aarseth, who made the NBODY4 
program for GRAPE-4 system available. I'm grateful 
to Makoto Taiji and Toshikazu Ebisuzaki, who de- 
veloped the GRAPE-4 system in collaboration with 
myself and Daiichiro Sugimoto. I thank Daiichiro 
Sugimoto for many stimulating discussions and help- 
ful comments on the manuscript. I also thank Pict 
Hut, Steve McMillan, Yoko Funato and Toshiyuki 
Fukushige for stimulating discussions. This work was 
supported by Grant-in-Aid for Specially Promoted 
Research (04102002) of the Ministry of Education, 
Science, Sports and Culture, Japan. 
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